home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / HENSA / MATHS / PLPLOT / PLPLOT.ZIP / src / plcont.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-06-30  |  22.1 KB  |  954 lines
  1. /* $Id: plcont.c,v 1.12 1994/06/30 18:22:02 mjl Exp $
  2.  * $Log: plcont.c,v $
  3.  * Revision 1.12  1994/06/30  18:22:02  mjl
  4.  * All core source files: made another pass to eliminate warnings when using
  5.  * gcc -Wall.  Lots of cleaning up: got rid of includes of math.h or string.h
  6.  * (now included by plplot.h), and other minor changes.  Now each file has
  7.  * global access to the plstream pointer via extern; many accessor functions
  8.  * eliminated as a result.
  9.  *
  10.  * Revision 1.11  1994/03/23  07:56:35  mjl
  11.  * Changed name of plcontf() to plfcont(), as part of the new 2d function
  12.  * plotter API.  plcontf is now a macro for backward compatibility.
  13.  *
  14.  * Replaced call to plexit() on simple (recoverable) errors with simply
  15.  * printing the error message (via plabort()) and returning.  Should help
  16.  * avoid loss of computer time in some critical circumstances (during a long
  17.  * batch run, for example).
  18. */
  19.  
  20. /*    plcont.c
  21.  
  22.     Contour plotter.
  23. */
  24.  
  25. #include "plplotP.h"
  26.  
  27. #ifdef MSDOS
  28. #pragma message("Microsoft programmers are sissies.")
  29. #pragma optimize("",off)
  30. #endif
  31.  
  32. /* Static function prototypes. */
  33.  
  34. static void
  35. plcntr(PLFLT (*plf2eval) (PLINT, PLINT, PLPointer),
  36.        PLPointer plf2eval_data,
  37.        PLINT nx, PLINT ny, PLINT kx, PLINT lx,
  38.        PLINT ky, PLINT ly, PLFLT flev, PLINT *iscan,
  39.        PLINT *ixstor, PLINT *iystor, PLINT nstor,
  40.        void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
  41.        PLPointer pltr_data);
  42.  
  43. static void
  44. pldrawcn(PLFLT (*plf2eval) (PLINT, PLINT, PLPointer),
  45.      PLPointer plf2eval_data,
  46.      PLINT nx, PLINT ny, PLINT kx, PLINT lx,
  47.      PLINT ky, PLINT ly, PLFLT flev, PLINT kcol, PLINT krow,
  48.      PLINT *p_kscan, PLINT *p_kstor, PLINT *iscan,
  49.      PLINT *ixstor, PLINT *iystor, PLINT nstor,
  50.      void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
  51.      PLPointer pltr_data);
  52.  
  53. static void
  54. plccal(PLFLT (*plf2eval) (PLINT, PLINT, PLPointer),
  55.        PLPointer plf2eval_data,
  56.        PLFLT flev, PLINT ix, PLINT iy,
  57.        PLINT ixg, PLINT iyg, PLFLT *dist);
  58.  
  59. static void 
  60. plr45 (PLINT *ix, PLINT *iy, PLINT isens);
  61.  
  62. static void 
  63. plr135 (PLINT *ix, PLINT *iy, PLINT isens);
  64.  
  65. /* Error flag for aborts */
  66.  
  67. static int error;
  68.  
  69. /*----------------------------------------------------------------------*\
  70.  * plf2eval2()
  71.  *
  72.  * Does a lookup from a 2d function array.  Array is of type (PLFLT **),
  73.  * and is column dominant (normal C ordering).
  74. \*----------------------------------------------------------------------*/
  75.  
  76. PLFLT
  77. plf2eval2(PLINT ix, PLINT iy, PLPointer plf2eval_data)
  78. {
  79.     PLFLT value;
  80.     PLfGrid2 *grid = (PLfGrid2 *) plf2eval_data;
  81.  
  82.     value = grid->f[ix][iy];
  83.  
  84.     return value;
  85. }
  86.  
  87. /*----------------------------------------------------------------------*\
  88.  * plf2eval()
  89.  *
  90.  * Does a lookup from a 2d function array.  Array is of type (PLFLT *),
  91.  * and is column dominant (normal C ordering).  You MUST fill the ny
  92.  * maximum array index entry in the PLfGrid struct.
  93. \*----------------------------------------------------------------------*/
  94.  
  95. PLFLT
  96. plf2eval(PLINT ix, PLINT iy, PLPointer plf2eval_data)
  97. {
  98.     PLFLT value;
  99.     PLfGrid *grid = (PLfGrid *) plf2eval_data;
  100.  
  101.     value = grid->f[ix * grid->ny + iy];
  102.  
  103.     return value;
  104. }
  105.  
  106. /*----------------------------------------------------------------------*\
  107.  * plf2evalr()
  108.  *
  109.  * Does a lookup from a 2d function array.  Array is of type (PLFLT *),
  110.  * and is row dominant (Fortran ordering).  You MUST fill the nx maximum
  111.  * array index entry in the PLfGrid struct.
  112. \*----------------------------------------------------------------------*/
  113.  
  114. PLFLT
  115. plf2evalr(PLINT ix, PLINT iy, PLPointer plf2eval_data)
  116. {
  117.     PLFLT value;
  118.     PLfGrid *grid = (PLfGrid *) plf2eval_data;
  119.  
  120.     value = grid->f[ix + iy * grid->nx];
  121.  
  122.     return value;
  123. }
  124.  
  125. /*----------------------------------------------------------------------*\
  126.  * void plcont()
  127.  *
  128.  * Draws a contour plot from data in f(nx,ny).  Is just a front-end to
  129.  * plfcont, with a particular choice for f2eval and f2eval_data.
  130. \*----------------------------------------------------------------------*/
  131.  
  132. void
  133. c_plcont(PLFLT **f, PLINT nx, PLINT ny, PLINT kx, PLINT lx,
  134.      PLINT ky, PLINT ly, PLFLT *clevel, PLINT nlevel,
  135.      void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
  136.      PLPointer pltr_data)
  137. {
  138.     PLfGrid2 grid;
  139.  
  140.     grid.f = f;
  141.     plfcont(plf2eval2, (PLPointer) &grid,
  142.         nx, ny, kx, lx, ky, ly, clevel, nlevel,
  143.         pltr, pltr_data);
  144. }
  145.  
  146. /*----------------------------------------------------------------------*\
  147.  * void plfcont()
  148.  *
  149.  * Draws a contour plot using the function evaluator f2eval and data
  150.  * stored by way of the f2eval_data pointer.  This allows arbitrary
  151.  * organizations of 2d array data to be used.
  152.  *
  153.  * The subrange of indices used for contouring is kx to lx in the x
  154.  * direction and from ky to ly in the y direction. The array of contour
  155.  * levels is clevel(nlevel), and "pltr" is the name of a function which
  156.  * transforms array indicies into world coordinates.
  157.  *
  158.  * Note that the fortran-like minimum and maximum indices (kx, lx, ky, ly)
  159.  * are translated into more C-like ones.  I've only kept them as they are
  160.  * for the plcontf() argument list because of backward compatibility.
  161. \*----------------------------------------------------------------------*/
  162.  
  163. void
  164. plfcont(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
  165.     PLPointer f2eval_data,
  166.     PLINT nx, PLINT ny, PLINT kx, PLINT lx,
  167.     PLINT ky, PLINT ly, PLFLT *clevel, PLINT nlevel,
  168.     void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
  169.     PLPointer pltr_data)
  170. {
  171.     PLINT i, mx, my, nstor, *heapc;
  172.  
  173.     mx = lx - kx + 1;
  174.     my = ly - ky + 1;
  175.  
  176.     if (kx < 1 || kx >= lx) {
  177.     plabort("plfcont: indices must satisfy  1 <= kx <= lx <= nx");
  178.     return;
  179.     }
  180.     if (ky < 1 || ky >= ly) {
  181.     plabort("plfcont: indices must satisfy  1 <= ky <= ly <= ny");
  182.     return;
  183.     }
  184.  
  185.     nstor = mx * my;
  186.     heapc = (PLINT *) malloc((size_t) (mx + 2 * nstor) * sizeof(PLINT));
  187.     if (heapc == NULL) {
  188.     plabort("plfcont: out of memory in heap allocation");
  189.     return;
  190.     }
  191.  
  192.     for (i = 0; i < nlevel; i++) {
  193.     plcntr(f2eval, f2eval_data,
  194.            nx, ny, kx-1, lx-1, ky-1, ly-1, clevel[i], &heapc[0],
  195.            &heapc[nx], &heapc[nx + nstor], nstor, pltr, pltr_data);
  196.  
  197.     if (error) {
  198.         error = 0;
  199.         goto done;
  200.     }
  201.     }
  202.  
  203.   done:
  204.     free((void *) heapc);
  205. }
  206.  
  207. /*----------------------------------------------------------------------*\
  208.  * void plcntr()
  209.  *
  210.  * The contour for a given level is drawn here.  Note iscan has nx
  211.  * elements. ixstor and iystor each have nstor elements.
  212. \*----------------------------------------------------------------------*/
  213.  
  214. static void
  215. plcntr(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
  216.        PLPointer f2eval_data,
  217.        PLINT nx, PLINT ny, PLINT kx, PLINT lx,
  218.        PLINT ky, PLINT ly, PLFLT flev, PLINT *iscan,
  219.        PLINT *ixstor, PLINT *iystor, PLINT nstor,
  220.        void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
  221.        PLPointer pltr_data)
  222. {
  223.     PLINT kcol, krow, kstor, kscan, l, ixt, iyt, jstor, next;
  224.  
  225.     /* Initialize memory pointers */
  226.  
  227.     kstor = 0;
  228.     kscan = 0;
  229.  
  230.     for (krow = ky; krow <= ly; krow++) {
  231.     for (kcol = kx + 1; kcol <= lx; kcol++) {
  232.  
  233.         /* Follow and draw a contour */
  234.  
  235.         pldrawcn(f2eval, f2eval_data,
  236.              nx, ny, kx, lx, ky, ly, flev, kcol, krow,
  237.              &kscan, &kstor, iscan, ixstor, iystor, nstor,
  238.              pltr, pltr_data);
  239.  
  240.         if (error)
  241.         return;
  242.     }
  243.  
  244.     /* Search of row complete */
  245.     /* Set up memory of next row in iscan and edit ixstor and iystor */
  246.  
  247.     if (krow < ny-1) {
  248.         jstor = 0;
  249.         kscan = 0;
  250.         next = krow + 1;
  251.         for (l = 1; l <= kstor; l++) {
  252.         ixt = ixstor[l - 1];
  253.         iyt = iystor[l - 1];
  254.  
  255.         /* Memory of next row into iscan */
  256.  
  257.         if (iyt == next) {
  258.             kscan = kscan + 1;
  259.             iscan[kscan - 1] = ixt;
  260.         }
  261.  
  262.         /* Retain memory of rows to come, and forget rest */
  263.  
  264.         else if (iyt > next) {
  265.             jstor = jstor + 1;
  266.             ixstor[jstor - 1] = ixt;
  267.             iystor[jstor - 1] = iyt;
  268.         }
  269.         }
  270.         kstor = jstor;
  271.     }
  272.     }
  273. }
  274.  
  275. /*----------------------------------------------------------------------*\
  276.  * void pldrawcn()
  277.  *
  278.  * Follow and draw a contour.
  279. \*----------------------------------------------------------------------*/
  280.  
  281. static void
  282. pldrawcn(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
  283.      PLPointer f2eval_data,
  284.      PLINT nx, PLINT ny, PLINT kx, PLINT lx,
  285.      PLINT ky, PLINT ly, PLFLT flev, PLINT kcol, PLINT krow,
  286.      PLINT *p_kscan, PLINT *p_kstor, PLINT *iscan,
  287.      PLINT *ixstor, PLINT *iystor, PLINT nstor,
  288.      void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
  289.      PLPointer pltr_data)
  290. {
  291.     PLINT iwbeg, ixbeg, iybeg, izbeg;
  292.     PLINT iboun, iw, ix, iy, iz, ifirst, istep, ixgo, iygo;
  293.     PLINT l, ixg, iyg, ia, ib;
  294.  
  295.     PLFLT dist, dx, dy, xnew, ynew, fxl, fxr;
  296.     PLFLT xlas = 0., ylas = 0., tpx, tpy, xt, yt;
  297.     PLFLT f1, f2, f3, f4, fcheck;
  298.  
  299.     /* Check if a contour has been crossed */
  300.  
  301.     fxl = f2eval(kcol-1, krow, f2eval_data);
  302.     fxr = f2eval(kcol, krow, f2eval_data);
  303.  
  304.     if (fxl < flev && fxr >= flev) {
  305.     ixbeg = kcol - 1;
  306.     iwbeg = kcol;
  307.     }
  308.     else if (fxr < flev && fxl >= flev) {
  309.     ixbeg = kcol;
  310.     iwbeg = kcol - 1;
  311.     }
  312.     else
  313.     return;
  314.  
  315.     iybeg = krow;
  316.     izbeg = krow;
  317.  
  318.     /* A contour has been crossed. */
  319.     /* Check to see if it is a new one. */
  320.  
  321.     for (l = 0; l < *p_kscan; l++) {
  322.     if (ixbeg == iscan[l])
  323.         return;
  324.     }
  325.  
  326.     for (iboun = 1; iboun >= -1; iboun -= 2) {
  327.  
  328.     /* Set up starting point and initial search directions */
  329.  
  330.     ix = ixbeg;
  331.     iy = iybeg;
  332.     iw = iwbeg;
  333.     iz = izbeg;
  334.     ifirst = 1;
  335.     istep = 0;
  336.     ixgo = iw - ix;
  337.     iygo = iz - iy;
  338.  
  339.     for (;;) {
  340.         plccal(f2eval, f2eval_data,
  341.            flev, ix, iy, ixgo, iygo, &dist);
  342.  
  343.         dx = dist * ixgo;
  344.         dy = dist * iygo;
  345.         xnew = ix + dx;
  346.         ynew = iy + dy;
  347.  
  348.         /* Has a step occured in search? */
  349.  
  350.         if (istep != 0) {
  351.         if (ixgo * iygo == 0) {
  352.  
  353.             /* This was a diagonal step, so interpolate missed */
  354.             /* point, rotating 45 degrees to get it */
  355.  
  356.             ixg = ixgo;
  357.             iyg = iygo;
  358.             plr45(&ixg, &iyg, iboun);
  359.             ia = iw - ixg;
  360.             ib = iz - iyg;
  361.             plccal(f2eval, f2eval_data,
  362.                flev, ia, ib, ixg, iyg, &dist);
  363.  
  364.             (*pltr) (xlas, ylas, &tpx, &tpy, pltr_data);
  365.  
  366.             plP_drawor(tpx, tpy);
  367.             dx = dist * ixg;
  368.             dy = dist * iyg;
  369.             xlas = ia + dx;
  370.             ylas = ib + dy;
  371.         }
  372.         else {
  373.             if (dist > 0.5) {
  374.             xt = xlas;
  375.             xlas = xnew;
  376.             xnew = xt;
  377.             yt = ylas;
  378.             ylas = ynew;
  379.             ynew = yt;
  380.             }
  381.         }
  382.         }
  383.         if (ifirst != 1) {
  384.         (*pltr) (xlas, ylas, &tpx, &tpy, pltr_data);
  385.         plP_drawor(tpx, tpy);
  386.         }
  387.         else {
  388.         (*pltr) (xnew, ynew, &tpx, &tpy, pltr_data);
  389.         plP_movwor(tpx, tpy);
  390.         }
  391.         xlas = xnew;
  392.         ylas = ynew;
  393.  
  394.         /* Check if the contour is closed */
  395.  
  396.         if (ifirst != 1 &&
  397.         ix == ixbeg && iy == iybeg && iw == iwbeg && iz == izbeg) {
  398.         (*pltr) (xlas, ylas, &tpx, &tpy, pltr_data);
  399.         plP_drawor(tpx, tpy);
  400.         return;
  401.         }
  402.         ifirst = 0;
  403.  
  404.         /* Now the rotation */
  405.  
  406.         istep = 0;
  407.         plr45(&ixgo, &iygo, iboun);
  408.         iw = ix + ixgo;
  409.         iz = iy + iygo;
  410.  
  411.         /* Check if out of bounds */
  412.  
  413.         if (iw < kx || iw > lx || iz < ky || iz > ly)
  414.         break;
  415.  
  416.         /* Has contact been lost with the contour? */
  417.  
  418.         if (ixgo * iygo == 0)
  419.         fcheck = f2eval(iw, iz, f2eval_data);
  420.         else {
  421.         f1 = f2eval(ix, iy, f2eval_data);
  422.         f2 = f2eval(iw, iz, f2eval_data);
  423.         f3 = f2eval(ix, iz, f2eval_data);
  424.         f4 = f2eval(iw, iy, f2eval_data);
  425.  
  426.         fcheck = MAX(f2, (f1 + f2 + f3 + f4) / 4.);
  427.         }
  428.  
  429.         if (fcheck < flev) {
  430.  
  431.         /* Yes, lost contact => step to new center */
  432.  
  433.         istep = 1;
  434.         ix = iw;
  435.         iy = iz;
  436.         plr135(&ixgo, &iygo, iboun);
  437.         iw = ix + ixgo;
  438.         iz = iy + iygo;
  439.  
  440.         /* And do the contour memory */
  441.  
  442.         if (iy == krow) {
  443.             *p_kscan = *p_kscan + 1;
  444.             iscan[*p_kscan - 1] = ix;
  445.         }
  446.         else if (iy > krow) {
  447.             *p_kstor = *p_kstor + 1;
  448.             if (*p_kstor > nstor) {
  449.             plabort("plfcont: heap exhausted");
  450.             error = 1;
  451.             return;
  452.             }
  453.             ixstor[*p_kstor - 1] = ix;
  454.             iystor[*p_kstor - 1] = iy;
  455.         }
  456.         }
  457.     }
  458.     /* Reach here only if boundary encountered - Draw last bit */
  459.  
  460.     (*pltr) (xnew, ynew, &tpx, &tpy, pltr_data);
  461.     plP_drawor(tpx, tpy);
  462.     }
  463. }
  464.  
  465. /*----------------------------------------------------------------------*\
  466.  * void plccal()
  467.  *
  468.  * Function to interpolate the position of a contour which is known to be
  469.  * next to ix,iy in the direction ixg,iyg. The unscaled distance along
  470.  * ixg,iyg is returned as dist.
  471. \*----------------------------------------------------------------------*/
  472.  
  473. static void
  474. plccal(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
  475.        PLPointer f2eval_data,
  476.        PLFLT flev, PLINT ix, PLINT iy,
  477.        PLINT ixg, PLINT iyg, PLFLT *dist)
  478. {
  479.     PLINT ia, ib;
  480.     PLFLT dbot, dtop, fmid;
  481.     PLFLT fxy, fab, fay, fxb, flow;
  482.  
  483.     ia = ix + ixg;
  484.     ib = iy + iyg;
  485.     fxy = f2eval(ix, iy, f2eval_data);
  486.     fab = f2eval(ia, ib, f2eval_data);
  487.     fxb = f2eval(ix, ib, f2eval_data);
  488.     fay = f2eval(ia, iy, f2eval_data);
  489.  
  490.     if (ixg == 0 || iyg == 0) {
  491.     dtop = flev - fxy;
  492.     dbot = fab - fxy;
  493.     *dist = 0.0;
  494.     if (dbot != 0.0)
  495.         *dist = dtop / dbot;
  496.     }
  497.     else {
  498.     fmid = (fxy + fab + fxb + fay) / 4.0;
  499.     *dist = 0.5;
  500.  
  501.     if ((fxy - flev) * (fab - flev) <= 0.) {
  502.  
  503.         if (fmid >= flev) {
  504.         dtop = flev - fxy;
  505.         dbot = fmid - fxy;
  506.         if (dbot != 0.0)
  507.             *dist = 0.5 * dtop / dbot;
  508.         }
  509.         else {
  510.         dtop = flev - fab;
  511.         dbot = fmid - fab;
  512.         if (dbot != 0.0)
  513.             *dist = 1.0 - 0.5 * dtop / dbot;
  514.         }
  515.     }
  516.     else {
  517.         flow = (fxb + fay) / 2.0;
  518.         dtop = fab - flev;
  519.         dbot = fab + fxy - 2.0 * flow;
  520.         if (dbot != 0.0)
  521.         *dist = 1. - dtop / dbot;
  522.     }
  523.     }
  524.     if (*dist > 1.)
  525.     *dist = 1.;
  526. }
  527.  
  528. /*----------------------------------------------------------------------*\
  529.  * Rotators
  530. \*----------------------------------------------------------------------*/
  531.  
  532. static void 
  533. plr45 (PLINT *ix, PLINT *iy, PLINT isens)
  534. {
  535.     PLINT ixx, iyy;
  536.  
  537.     ixx = *ix - isens * (*iy);
  538.     iyy = *ix * isens + *iy;
  539.     *ix = ixx / MAX(1, ABS(ixx));
  540.     *iy = iyy / MAX(1, ABS(iyy));
  541. }
  542.  
  543. static void 
  544. plr135 (PLINT *ix, PLINT *iy, PLINT isens)
  545. {
  546.     *ix = -*ix;
  547.     *iy = -*iy;
  548.     plr45(ix, iy, isens);
  549. }
  550.  
  551. /*----------------------------------------------------------------------*\
  552.  * pltr0()
  553.  *
  554.  * Identity transformation.
  555. \*----------------------------------------------------------------------*/
  556.  
  557. void
  558. pltr0(PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data)
  559. {
  560.     *tx = x;
  561.     *ty = y;
  562. }
  563.  
  564. /*----------------------------------------------------------------------*\
  565.  * pltr1()
  566.  *
  567.  * Does linear interpolation from singly dimensioned coord arrays.
  568.  *
  569.  * Just abort for now if coordinates are out of bounds (don't think it's
  570.  * possible, but if so we could use linear extrapolation).
  571. \*----------------------------------------------------------------------*/
  572.  
  573. void
  574. pltr1(PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data)
  575. {
  576.     PLINT ul, ur, vl, vr;
  577.     PLFLT du, dv;
  578.     PLFLT xl, xr, yl, yr;
  579.  
  580.     PLcGrid *grid = (PLcGrid *) pltr_data;
  581.     PLFLT *xg = grid->xg;
  582.     PLFLT *yg = grid->yg;
  583.     PLINT nx = grid->nx;
  584.     PLINT ny = grid->ny;
  585.  
  586.     ul = (PLINT) x;
  587.     ur = ul + 1;
  588.     du = x - ul;
  589.  
  590.     vl = (PLINT) y;
  591.     vr = vl + 1;
  592.     dv = y - vl;
  593.  
  594.     if (x < 0 || x > nx - 1 || y < 0 || y > ny - 1) {
  595.     plexit("pltr1: Invalid coordinates");
  596.     }
  597.  
  598. /* Look up coordinates in row-dominant array.
  599.  * Have to handle right boundary specially -- if at the edge, we'd
  600.  * better not reference the out of bounds point. 
  601.  */
  602.  
  603.     xl = xg[ul];
  604.     yl = yg[vl];
  605.  
  606.     if (ur == nx) {
  607.     *tx = xl;
  608.     }
  609.     else {
  610.     xr = xg[ur];
  611.     *tx = xl * (1 - du) + xr * du;
  612.     }
  613.     if (vr == ny) {
  614.     *ty = yl;
  615.     }
  616.     else {
  617.     yr = yg[vr];
  618.     *ty = yl * (1 - dv) + yr * dv;
  619.     }
  620. }
  621.  
  622. /*----------------------------------------------------------------------*\
  623.  * pltr2()
  624.  *
  625.  * Does linear interpolation from doubly dimensioned coord arrays (column
  626.  * dominant, as per normal C 2d arrays).
  627.  *
  628.  * This routine includes lots of checks for out of bounds.  This would
  629.  * occur occasionally due to some bugs in the contour plotter (now fixed).
  630.  * If an out of bounds coordinate is obtained, the boundary value is
  631.  * provided along with a warning.  These checks should stay since no harm
  632.  * is done if if everything works correctly.
  633. \*----------------------------------------------------------------------*/
  634.  
  635. void
  636. pltr2(PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data)
  637. {
  638.     PLINT ul, ur, vl, vr;
  639.     PLFLT du, dv;
  640.     PLFLT xll, xlr, xrl, xrr;
  641.     PLFLT yll, ylr, yrl, yrr;
  642.     PLFLT xmin, xmax, ymin, ymax;
  643.  
  644.     PLcGrid2 *grid = (PLcGrid2 *) pltr_data;
  645.     PLFLT **xg = grid->xg;
  646.     PLFLT **yg = grid->yg;
  647.     PLINT nx = grid->nx;
  648.     PLINT ny = grid->ny;
  649.  
  650.     ul = (PLINT) x;
  651.     ur = ul + 1;
  652.     du = x - ul;
  653.  
  654.     vl = (PLINT) y;
  655.     vr = vl + 1;
  656.     dv = y - vl;
  657.  
  658.     xmin = 0;
  659.     xmax = nx - 1;
  660.     ymin = 0;
  661.     ymax = ny - 1;
  662.  
  663.     if (x < xmin || x > xmax || y < ymin || y > ymax) {
  664.     plwarn("pltr2: Invalid coordinates");
  665.     if (x < xmin) {
  666.  
  667.         if (y < ymin) {
  668.         *tx = xg[0][0];
  669.         *ty = yg[0][0];
  670.         }
  671.         else if (y > ymax) {
  672.         *tx = xg[0][ny-1];
  673.         *ty = yg[0][ny-1];
  674.         }
  675.         else {
  676.         xll = xg[0][vl];
  677.         yll = yg[0][vl];
  678.         xlr = xg[0][vr];
  679.         ylr = yg[0][vr];
  680.  
  681.         *tx = xll * (1 - dv) + xlr * (dv);
  682.         *ty = yll * (1 - dv) + ylr * (dv);
  683.         }
  684.     }
  685.     else if (x > xmax) {
  686.  
  687.         if (y < ymin) {
  688.         *tx = xg[nx-1][0];
  689.         *ty = yg[nx-1][0];
  690.         }
  691.         else if (y > ymax) {
  692.         *tx = xg[nx-1][ny-1];
  693.         *ty = yg[nx-1][ny-1];
  694.         }
  695.         else {
  696.         xll = xg[nx-1][vl];
  697.         yll = yg[nx-1][vl];
  698.         xlr = xg[nx-1][vr];
  699.         ylr = yg[nx-1][vr];
  700.  
  701.         *tx = xll * (1 - dv) + xlr * (dv);
  702.         *ty = yll * (1 - dv) + ylr * (dv);
  703.         }
  704.     }
  705.     else {
  706.         if (y < ymin) {
  707.         xll = xg[ul][0];
  708.         xrl = xg[ur][0];
  709.         yll = yg[ul][0];
  710.         yrl = yg[ur][0];
  711.  
  712.         *tx = xll * (1 - du) + xrl * (du);
  713.         *ty = yll * (1 - du) + yrl * (du);
  714.         }
  715.         else if (y > ymax) {
  716.         xlr = xg[ul][ny-1];
  717.         xrr = xg[ur][ny-1];
  718.         ylr = yg[ul][ny-1];
  719.         yrr = yg[ur][ny-1];
  720.  
  721.         *tx = xlr * (1 - du) + xrr * (du);
  722.         *ty = ylr * (1 - du) + yrr * (du);
  723.         }
  724.     }
  725.     }
  726.  
  727. /* Normal case.
  728.  * Look up coordinates in row-dominant array.
  729.  * Have to handle right boundary specially -- if at the edge, we'd
  730.  * better not reference the out of bounds point. 
  731.  */
  732.  
  733.     else {
  734.  
  735.     xll = xg[ul][vl];
  736.     yll = yg[ul][vl];
  737.  
  738. /* ur is out of bounds */
  739.  
  740.     if (ur == nx && vr < ny) {
  741.  
  742.         xlr = xg[ul][vr];
  743.         ylr = yg[ul][vr];
  744.  
  745.         *tx = xll * (1 - dv) + xlr * (dv);
  746.         *ty = yll * (1 - dv) + ylr * (dv);
  747.     }
  748.  
  749. /* vr is out of bounds */
  750.  
  751.     else if (ur < nx && vr == ny) {
  752.  
  753.         xrl = xg[ur][vl];
  754.         yrl = yg[ur][vl];
  755.  
  756.         *tx = xll * (1 - du) + xrl * (du);
  757.         *ty = yll * (1 - du) + yrl * (du);
  758.     }
  759.  
  760. /* both ur and vr are out of bounds */
  761.  
  762.     else if (ur == nx && vr == ny) {
  763.  
  764.         *tx = xll;
  765.         *ty = yll;
  766.     }
  767.  
  768. /* everything in bounds */
  769.  
  770.     else {
  771.  
  772.         xrl = xg[ur][vl];
  773.         xlr = xg[ul][vr];
  774.         xrr = xg[ur][vr];
  775.  
  776.         yrl = yg[ur][vl];
  777.         ylr = yg[ul][vr];
  778.         yrr = yg[ur][vr];
  779.  
  780.         *tx = xll * (1 - du) * (1 - dv) + xlr * (1 - du) * (dv) +
  781.         xrl * (du) * (1 - dv) + xrr * (du) * (dv);
  782.  
  783.         *ty = yll * (1 - du) * (1 - dv) + ylr * (1 - du) * (dv) +
  784.         yrl * (du) * (1 - dv) + yrr * (du) * (dv);
  785.     }
  786.     }
  787. }
  788.  
  789. /*----------------------------------------------------------------------*\
  790.  * pltr2p()
  791.  *
  792.  * Just like pltr2() but uses pointer arithmetic to get coordinates from
  793.  * 2d grid tables.  This form of grid tables is compatible with those from
  794.  * PLPLOT 4.0.  The grid data must be pointed to by a PLcGrid structure.
  795. \*----------------------------------------------------------------------*/
  796.  
  797. void
  798. pltr2p(PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data)
  799. {
  800.     PLINT ul, ur, vl, vr;
  801.     PLFLT du, dv;
  802.     PLFLT xll, xlr, xrl, xrr;
  803.     PLFLT yll, ylr, yrl, yrr;
  804.     PLFLT xmin, xmax, ymin, ymax;
  805.  
  806.     PLcGrid *grid = (PLcGrid *) pltr_data;
  807.     PLFLT *xg = grid->xg;
  808.     PLFLT *yg = grid->yg;
  809.     PLINT nx = grid->nx;
  810.     PLINT ny = grid->ny;
  811.  
  812.     ul = (PLINT) x;
  813.     ur = ul + 1;
  814.     du = x - ul;
  815.  
  816.     vl = (PLINT) y;
  817.     vr = vl + 1;
  818.     dv = y - vl;
  819.  
  820.     xmin = 0;
  821.     xmax = nx - 1;
  822.     ymin = 0;
  823.     ymax = ny - 1;
  824.  
  825.     if (x < xmin || x > xmax || y < ymin || y > ymax) {
  826.     plwarn("pltr2p: Invalid coordinates");
  827.     if (x < xmin) {
  828.  
  829.         if (y < ymin) {
  830.         *tx = *xg;
  831.         *ty = *yg;
  832.         }
  833.         else if (y > ymax) {
  834.         *tx = *(xg + (ny - 1));
  835.         *ty = *(yg + (ny - 1));
  836.         }
  837.         else {
  838.         ul = 0;
  839.         xll = *(xg + ul * ny + vl);
  840.         yll = *(yg + ul * ny + vl);
  841.         xlr = *(xg + ul * ny + vr);
  842.         ylr = *(yg + ul * ny + vr);
  843.  
  844.         *tx = xll * (1 - dv) + xlr * (dv);
  845.         *ty = yll * (1 - dv) + ylr * (dv);
  846.         }
  847.     }
  848.     else if (x > xmax) {
  849.  
  850.         if (y < ymin) {
  851.         *tx = *(xg + (ny - 1) * nx);
  852.         *ty = *(yg + (ny - 1) * nx);
  853.         }
  854.         else if (y > ymax) {
  855.         *tx = *(xg + (ny - 1) + (nx - 1) * ny);
  856.         *ty = *(yg + (ny - 1) + (nx - 1) * ny);
  857.         }
  858.         else {
  859.         ul = nx - 1;
  860.         xll = *(xg + ul * ny + vl);
  861.         yll = *(yg + ul * ny + vl);
  862.         xlr = *(xg + ul * ny + vr);
  863.         ylr = *(yg + ul * ny + vr);
  864.  
  865.         *tx = xll * (1 - dv) + xlr * (dv);
  866.         *ty = yll * (1 - dv) + ylr * (dv);
  867.         }
  868.     }
  869.     else {
  870.         if (y < ymin) {
  871.         vl = 0;
  872.         xll = *(xg + ul * ny + vl);
  873.         xrl = *(xg + ur * ny + vl);
  874.         yll = *(yg + ul * ny + vl);
  875.         yrl = *(yg + ur * ny + vl);
  876.  
  877.         *tx = xll * (1 - du) + xrl * (du);
  878.         *ty = yll * (1 - du) + yrl * (du);
  879.         }
  880.         else if (y > ymax) {
  881.         vr = ny - 1;
  882.         xlr = *(xg + ul * ny + vr);
  883.         xrr = *(xg + ur * ny + vr);
  884.         ylr = *(yg + ul * ny + vr);
  885.         yrr = *(yg + ur * ny + vr);
  886.  
  887.         *tx = xlr * (1 - du) + xrr * (du);
  888.         *ty = ylr * (1 - du) + yrr * (du);
  889.         }
  890.     }
  891.     }
  892.  
  893. /* Normal case.
  894.  * Look up coordinates in row-dominant array.
  895.  * Have to handle right boundary specially -- if at the edge, we'd
  896.  * better not reference the out of bounds point. 
  897.  */
  898.  
  899.     else {
  900.  
  901.     xll = *(xg + ul * ny + vl);
  902.     yll = *(yg + ul * ny + vl);
  903.  
  904. /* ur is out of bounds */
  905.  
  906.     if (ur == nx && vr < ny) {
  907.  
  908.         xlr = *(xg + ul * ny + vr);
  909.         ylr = *(yg + ul * ny + vr);
  910.  
  911.         *tx = xll * (1 - dv) + xlr * (dv);
  912.         *ty = yll * (1 - dv) + ylr * (dv);
  913.     }
  914.  
  915. /* vr is out of bounds */
  916.  
  917.     else if (ur < nx && vr == ny) {
  918.  
  919.         xrl = *(xg + ur * ny + vl);
  920.         yrl = *(yg + ur * ny + vl);
  921.  
  922.         *tx = xll * (1 - du) + xrl * (du);
  923.         *ty = yll * (1 - du) + yrl * (du);
  924.     }
  925.  
  926. /* both ur and vr are out of bounds */
  927.  
  928.     else if (ur == nx && vr == ny) {
  929.  
  930.         *tx = xll;
  931.         *ty = yll;
  932.     }
  933.  
  934. /* everything in bounds */
  935.  
  936.     else {
  937.  
  938.         xrl = *(xg + ur * ny + vl);
  939.         xlr = *(xg + ul * ny + vr);
  940.         xrr = *(xg + ur * ny + vr);
  941.  
  942.         yrl = *(yg + ur * ny + vl);
  943.         ylr = *(yg + ul * ny + vr);
  944.         yrr = *(yg + ur * ny + vr);
  945.  
  946.         *tx = xll * (1 - du) * (1 - dv) + xlr * (1 - du) * (dv) +
  947.         xrl * (du) * (1 - dv) + xrr * (du) * (dv);
  948.  
  949.         *ty = yll * (1 - du) * (1 - dv) + ylr * (1 - du) * (dv) +
  950.         yrl * (du) * (1 - dv) + yrr * (du) * (dv);
  951.     }
  952.     }
  953. }
  954.